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Two first-principles simulation techniques, path integral Monte Carlo (PIMC) and density func- 
tional molecular dynamics (DFT-MD), are applied to study hot, dense helium in the density- 
temperature range of 0.387 - 5.35 gem"'' and 500 K - 1.28x10* K. One coherent equation of 
state (EOS) is derived by combining DFT-MD data at lower temperatures with PIMC results at 
higher temperatures. Good agreement between both techniques is found in an intermediate temper- 
ature range. For the highest temperatures, the PIMC results converge to the Debye-Hiickel limiting 
law. In order derive the entropy, a thermodynamically consistent free energy fit is introduced that 
reproduces the internal energies and pressure derived from the first-principles simulations. The 
equation of state is presented in form of a table as well as a fit and is compared with chemical 
models. In addition, the structure of the fiuid is analyzed using pair correlation functions. Shock 
Hugoniot curves are compared with recent laser shock wave experiments. 
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I. INTRODUCTION 

After hydrogen, helium is the second most common cl- 
ement in the universe. While it rarely occurs in pure form 
in nature, it is an endmember of hydrogen-helium mix- 
tures (HHM) that are the prevalent component in solar 
and extrasolar giant gas planets. The characterization of 
helium's properties at extreme temperature and pressure 
conditions is therefore important to study planetary in- 
teriors and especially relevant for answering the question 
whether HHM phase-separate in giant planet intcriors^'^. 
In most planetary models, the equation of state (EOS) of 
HHM was inferred from the linear mixing approximation 
at constant pressure and temperature using the EOSs of 
pure hydrogen and helium. The latter is the central topic 
of this article. 

Hydrogen and helium share some common properties. 
Both are very light and exhibit rich quantum proper- 
ties at low temperature. More importantly for this pa- 
per, the helium atom and the deuterium molecule have 
similar masses and both have two elemental excitation 
mechanisms that determine their behavior at high tem- 
perature. The helium atom has two ionization stages 
while deuterium molecules can dissociate and the result- 
ing atoms can be ionized. However, helium is without 
question simpler to characterize at high pressure. The 
crystal structure is hexagonal closed-packed under most 
(P, T) conditions^'"' while in solid hydrogen, different de- 
grees of molecular rotational ordering lead to several 
phases that deviate from the h.c.p. structure. Hydro- 
gen is expected to turn metallic at a few hundred GPa 
while a much larger bandgap must be closed in helium, 
which is predicted to occur above 10 000 eFax's. 

Given the relative simplicity of helium's high pressure 
properties one expects that there would be less of a con- 
troversy in the EOS than for hydrogen. This makes he- 
lium a good material to test novel experimental and the- 
oretical approaches. For hydrogen, the first laser shock 
experiments that reached megabar pressures had pre- 



dicted that the material would be highly compressible un- 
der shock conditions and reach densities six times higher 
than initial state^i^. Later experiments2.ii£iiiii- showed 
reduced compression ratios close to 4.3, which were in 
good agreement with first-principles calculation !] "'^i^^'^^ . 
In the case of helium, there is very good agreement be- 
tween the early shock experiments by Nellis et al.^^ and 
first-principles calculation^^. 

Recently the first laser shock experiments were per- 
formed on precompressed helium samples^^. The mea- 
surements confirmed the theoretically predicted trendi^ 
that the shock compression ratio is reduced with increas- 
ing precompression. However, there is discrepancy in the 
magnitude of the compression. Shock measurementsiS 
without precompression showed compression ratios of 
about 6 while first-principles simulation^'' predicted only 
5.24(4). The discrepancy between theoretical and experi- 
mental predictions is reduced for higher precompressions. 
For samples that were precompressed to 3.4-fold the am- 
bient density, theory and experiment are in agreement. 

The properties of fluid helium change from hard-sphere 
liquid at low pressure and temperature to ultimately a 
two-component plasma, after ionization has occurred at 
high pressure and temperature. The associated insulator- 
to-metal transition has been the topic of three recent the- 
oretical studies that all relied on DFT methods. Kietz- 
mann et ali^ studied the rise in electrical conductivity 
using the Kubo-Greenwood formula and compared with 
shock-wave experiments by Ternovoi et ali^. Kowalski 
et al.^ studied dense helium in order to characterize the 
atmosphere of white dwarfs. This paper went beyond 
the generalized gradient approximation by considering 
hybrid functionals. Stixrude and Jeanloz'^i studied the 
band gap closure in the dense fluid helium at over a wide 
range of densities including conditions of giant planet in- 
teriors. 

This article provides the EOS for fluid helium over a 
wide range of temperatures (500K-1.28xlO^K) and den- 
sities (0.387-5.35 gcm"'^ corresponding to a Wigner-Seitz 
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FIG. 1: The upper panel shows the finite size dependence of 
the pressure as function of the number atoms, A^, as pre- 
dicted from PIMC simulations with free-particle nodes at 
T=125 000K and Vg — 1.86. The lower panel compares the 
finite size dependence of DFT-MD simulations and classi- 
cal Monte Carlo calculations using the Aziz pair potential 
at T^IOOOOK and = 2.4. 



II. METHODS 

Path integral Monte Carlc^^ is the most appropri- 
ate and efficient first-principles simulation techniques for 
quantum system with thermal excitations. Electrons and 
nuclei are treated equally as paths, although the zero- 
point motion of the nuclei as well as exchange effects are 
negligible for the temperatures under consideration. The 
Coulomb interaction between electrons and nuclei is in- 
troduced using pair density matrices that we derived us- 
ing the eigenstates of the two-body Coulomb probleni^^. 
The periodic images were treated using an optimized 
Ewald break-up— that we applied to the pair action^S,. 
The explicit treatment of electrons as paths leads to the 
fermion sign problem, which requires one to introduce 
the only uncontrolled approximation in this method, the 
fixed node approximation2Si21. We use the nodes from 
the free-particle density matrix and from a variational 
density matrix^^. Besides this approximation, all corre- 
lation effects are included in PIMC, which for example 
leads an exact treatment of the isolated helium atom. 

The DFT-MD simulations were performed with ei- 
ther the CPMD code^^ using local TrouUier-Martins 
norm-conserving pseudopotentials^ or with the Vi- 
enna ab initio simulation package^! using the projec- 
tor augmented-wave method^. The nuclei were propa- 
gated using Born-Oppenheimer molecular dynamics with 
forces derived from either the electronic ground state 
or by including thermally excited electronic states when 
needed. Exchange-correlation effects were described 
by the Perdew-Burke-Ernzerhof generalized gradient ap- 
proximation'^^. The electronic wavefunctions were ex- 
panded in a plane- wave basis with energy cut-off of 30-50 
Hartrees. Most simulations were performed with iV=64 
using F point sampling of the Brillioun zone. An analysis 
of finite size effect is presented in the following section. 



radius interval of rs=2. 4-1.0 where ^nr^ ~ V/N^) by 
combining two first-principles simulation methods, path 
integral Monte Carlo (PIMC) at higher temperatures 
with density functional molecular dynamics (DFT-MD) 
at lower temperatures. The temperature range was sig- 
nificantly extended compared to our earlier worki^ that 
focused exclusively on shock properties alone. Here, the 
region of validity of both first-principles methods is ana- 
lyzed and good agreement for EOS at intermediate tem- 
peratures is demonstrated. The PIMC calculations have 
been extended to much higher temperatures until good 
agreement with the Debye-Hiickel limiting law is found. 
In the density interval under consideration, the entire 
EOS of nonrelativistic, fluid helium has been mapped 
out from first principles. All EOS data are combined into 
one thermodynamically consistent fit for the free energy 
and the entropy is derived. The structure of the fluid is 
analyzed using pair correlation functions and, finally, ad- 
ditional results for shock Hugoniot curves are presented. 



III. EQUATION OF STATE 

An analysis of finite size dependence of the EOS results 
is important since all simulation are perform with a fi- 
nite number of particle in periodic boundary conditions. 
Figure [1] gives two examples for the finite size analysis 
that we have performed at various (T, p) conditions. At 
10 000 K and = 2.4, helium can be characterized as a 
hard-sphere fiuid. The artificial periodicity of the nuclei 
dominate the finite size errors. Simulations with iV = 64 
atoms are sufficiently accurate for the purpose of this 
study. The DFT-MD results also agree surprisingly well 
with classical Monte Carlo calculation using the Aziz pair 
potential"^'*, which explains why both methods give fairly 
similar Hugoniot curves as long as electronic excitations 
are not importanl^. 

The upper panel of Fig. [T] shows PIMC results for 
125 000 K where a substantial part of the pressure comes 
from excited electrons. They are still coupled to the mo- 
tion of the nuclei, which leads to effective screening. In 
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FIG. 2: Comparison of the relative excess pressure de- 
rived from PIMC (solid lines) and DFT-MD. The dashed and 
the dash-dotted lines show results DFT-MD simulation with 
and without the consideration of thermally excited electronic 
states, respectively. 



consequence, the finite size dependence of the pressure 
is reduced significantly, and a simulation with = 16 
atoms exhibit a finite size error of only 1% compared with 
3% at lower temperature. This is reason why PIMC simu- 
lations with 16 atoms already gave a fairly accurate shock 
Hugoniot curveii. However, most PIMC results reported 
in Tab. U were obtained with 32 atoms and some with 57 
atoms. Already 32 atoms lead to well converged pressures 
unless one is interested at very high temperature above 
lO^K where all atoms are ionized and the coupling is 
very weak. Although the total pressure is dominated by 
the kinetic term, the excess pressure shows an increased 
finite size dependence that requires simulation with 57 
atoms in some cases. In general, the weak-coupling limit 
is difficult to study with finite-size simulations'^^. Also 
at very high density beyond the range considered here, 
electrons approach the limit of an ideal Fermi gas and 
form a rigid background. The remaining Coulombic sub- 
system of ions is expected to require simulations with 
several hundreds of particles^. In this regard, the elec- 
tronic screening makes our simulations affordable. 

Figure O compares the pressures obtained from PIMC 
and DFT-MD simulations for several density. Above 
20 000 K, excited electronic state become important. 
Both first-principles method are in very good agreement, 
which is foundation for the coherent EOS reported in this 
paper. Reasonably good agreement between PIMC and 
DFT-MD was reported for hydrogen earlier—. Figure [5] 
is a stringent test because it compares only the pressure 
contributions that result from the particle interactions. 
When one removes the ideal gas contributions, Pq, has 
one a bit of a choice for the corresponding noninteracting 
system. At very high temperature, one wants to com- 
pare with an ideal Fermi gas of electrons and nuclei. At 
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FIG. 3: Excess internal energy per electron relative to the 
ideal plasma model at a density of rs = 1.86. The circles show 
PIMC results. In the upper panel, the open squares and dia- 
monds show DFT-MD results with and without thermal pop- 
ulation of excited electronic states; respectively. The filled 
squares shows DFT-MD results corrected by constant shift 
corresponding to the DFT error of the isolated helium atom. 
In lower panel, PIMC results are compared with the Debye 
model. 



low temperature, however, one prefers comparing with an 
gas of noninteracting atoms. To combine these to lim- 
iting cases, we construct a simple chemical model that 
includes neutral atoms. He"*" and He^^ ions as well as 
free electrons. The ionization state is derived from the 
Saha equilibrium using the following binding energies, 
Ehs = -79.0 eV, EHe+ = -54.4 eV. Besides the binding 
energies, no other interaction are considered. 

This approach smoothly connects the low- and high- 
temperature limits. However for the correct interpreta- 
tion of the presented graphs, it should be pointed out 
that the pressures and energies of the ideal model de- 
pend on the Saha ionization equilibrium. If the ideal 
system exhibits a higher degree of ionization than the 
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TABLE L EOS table with pressures, internal and free ener- 
gies per electron derived from (a) DFT-MD with 64 atoms (a 
uniform AE/N^ = —0.4909 Ha correction was added to ac- 
count for missing DFT correlation energy in the helium atom), 
PIMC with (b) 32 atoms, (c) PIMC with 57 atoms, and (d) 
Debye-Hiickel limiting law. The numbers in brackets indi- 
cate the statistical uncertainties of the DFT-MD and PIMC 
simulations for the corresponding number trailing digits. 





T(K) 


P (GPa) 


tp / AT f'tJn\ 


t /l\e (tia.) 




Z.4 


500 


1.420(10) 


1 /1/10Q7'Q/'7\ 
-i.44yo io[i) 


-i.4004 




1000 


2.045(14) 




-1.40100 


2.4° 


3000 


4.69(3) 




-1.49126 


2.4" 


5000 


6.98(4) 


-1.43727(3) 


-1.52534 


2.4" 


10000 


12.49(4) 


-1.42395(5) 


-1.61873 


2.4" 


20000 


22.19(8) 


-1.39427(12) 


-1.82386 


2.4" 


40000 


43.37(11) 


-1.2997(2) 


-2.28643 


2.4" 


60000 


68.27(10) 


-1.1748(2) 


-2.80627 


2.4" 


80000 


96.93(12) 


-1.02525(7) 


-3.37236 


2.4'' 


125000 


172.3(6) 


-0.6667(17) 


-4.77369 


2.4'' 


250000 


445.7(6) 


0.477(2) 


-9.31702 


2.4'' 


333333 


651.4(9) 


1.237(3) 


-12.6707 


2.4'' 


500000 1067.7(1.0) 


2.634(3) 


-19.922 


2.4'' 


571428 


1249.9(9) 


3.216(3) 


-23.1952 


2.4'' 


666667 


1484.2(5) 


3.9507(16) 


-27.6612 


2.4'' 


800000 1815.5(1.2) 


4.972(4) 


-34.0708 


2.4'' 


1x10** 


2308.4(7) 


6.470(2) 


-43.9954 


2.4'' 


2x10** 


4745.2(8) 


13.754(2) 


-97.4249 


2.4'' 


4x10*^ 


9587.6(1.2) 


28.102(4) 


-213.949 


2.4'* 


SxlO'' 


19253 


56.72 


-466.803 


2.4'* 


16x10^ 


38577 


113.80 


-1013.36 


2.4'* 


32x10^ 


77205 


227.87 


-2184.48 


2.4'* 


64x10*^ 


154445 


455.92 


-4683.43 


2.4'* 


128x10*^ 


308916 


911.97 


-10002.6 


2.4'* 


256x10*^ 


617849 


1824.04 


-21267.5 


2.4'* 


512x10^ 


1235711 


3648.14 


-45052.4 


2.4'* 


1024x10** 


2471430 


7296.32 


-95193.5 


2.4'* 


2048x10^ 


4942866 


14592.67 


-200489 



simulation results, then this alone can lead to negative 
excess pressures and energies, which one would normally 
attribute exclusively to the interaction of free electrons 
and ions. This effect becomes clear in Fig. [3] where even 
the DFT-MD results without excited electrons exhibit a 
negative excess internal energy at 80 000 K. 

Figure [3] exhibits the missing correlation energies in 
DFT GGA, which underestimates the binding energy in 
the atom by A£'o=0.98 eV. This is the main reason for 
the disagreement with the ideal plasma model in the low- 
temperature limit (1.20 eV/atom), while internal energy 
increase due to the compression to = 1.86 amounts 
only to 0.22 eV/atom. To correct for missing correlation 
energy, we applied a uniform correction of —AEq to all 
DFT results discussed later. We may assume that the 
correction to DFT depends only weakly on temperature 
and density. To determine its precise amount more ac- 
curately is difficult and goes beyond the scope of this 
article. 

Despite this DFT insufficiency, one finds reasonably 
good agreement in the internal energies reported by 



TABLE II; Table U continued. 



rs 


r(K) 


P (GPa) 


E/Ne (Ha) 


J~l/A7 /TT \ 

F/Ne (Ha) 


2.2" 


500 


2.74(2) 


-1.449495(13) 


-1.45451 


2.2" 


1000 


3.77(3) 


-1.44787(2) 


-1.4601 


2.2" 


3000 


7.59(3) 


-1.44186(2) 


-1.48855 


2.2" 


5000 


10.81(6) 


-1.43615(3) 


-1.5214 


2.2 
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9 9^ 
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2.2" 
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-2 2635 


2.2" 


60000 


92.10(11) 


-1.17937(16) 


-2.77064 


2.2" 


80000 


129.52(10) 


-1.03399(12) 


-3.32222 


2.2'' 


125000 


223.9(7) 


-0.6962(16) 


-4.68439 


2.2'' 


250000 


569.6(7) 


0.3971(16) 


-9.08923 


2.2'' 


500000 


1371.3(6) 


2.5403(14) 


-19.3777 


2.2'' 


IxlO*" 


2981.1(7) 


6.3896(17) 


-42.807 


2.2'' 


2x10^ 


6148.3(8) 


13.694(2) 


-95 


2.2'' 


4x10^ 


12438.6(1.7) 


28.059(4) 


-209.039 


2.2'* 


SxlO** 


24982 


56.67 


-456.924 


2.2'* 


16x10*" 


50075 


113.77 


-993.605 


2.2'* 


32x10** 


100227 


227.85 


-2144.93 


2.2'* 


64x10** 


200508 


455.91 


-4604.15 


2.2'* 
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1 fin/19Q7 
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-444ZU.0 


2.2'* 


1024 X 10 


oono c o T 

3208587 


7296.31 


nonoo A 

-93928.4 


2.2'* 


2048x10 


6417184 


14592.66 


-197960 


2" 


500 


6.101(13) 


-1.448584(6) 


-1.45297 


2" 


1000 


7.59(2) 


-1.446835(11) 


-1.45806 




3000 


13.57(5) 


-1.44022(3) 


-1.48466 


2" 


5000 


18.07(7) 


-1.43433(4) 


-1.51606 


I 


lUUUU 


28.02(iij 


1 A OA1 O ^ Q\ 

-i.42013(8j 


-1.60363 
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9nnnn 

zuuuu 






1 7Q771 


qo, 


4nnnn 


84 79n 9") 


-1 30039f 1 5"! 


-9 93673 


qo, 


finnnn 


1 98 66n 4") 


-1 18357n9'l 


-9 79Q86 


2" 


80000 


178.84(19) 


-1.0433(2) 


-3.26574 


2' 


125000 


297.4(7) 


-0.7291(12) 


-4.58614 


2'' 


250000 


745.9(7) 


0.3112(13) 


-8.84403 


2* 


500000 


1800.9(8) 


2.4269(14) 


-18.7895 


2* 


1x10*^ 


3941.2(1.0) 


6.2898(17) 


-41.5111 


2'' 


2x10*^ 


8163.1(1.5) 


13.621(3) 


-92.3417 


2'' 


4x10** 


16544(2) 


28.009(4) 


-203.653 


2'* 


8x10** 


33228 


56.61 


-446.089 


2'* 


16x10** 


66633 


113.73 


-971.926 


2^ 


32x 10^ 


1 33390 


227.82 


-9101 53 


2'* 


64x10** 


266868 


455.88 


-4517.11 


2'* 


128x10** 


533797 


911.95 


-9669.7 


2'* 


256x10** 


1067636 


1824.02 


-20604.1 


2'* 


512x10** 


2135303 


3648.12 


-43727.2 


2'* 


1024x10** 


4270628 


7296.31 


-92540.1 


2'* 


2048x10** 


8541271 


14592.66 


-195186 



PIMC and DFT-MD. Figure [3] shows that both meth- 
ods report very similar increases in the energy resulting 
from thermal excitations of electrons, which is the basis 
for constructing one EOS table. 

In order to explore the agreement between PIMC and 
DFT-MD in more detail we resort to pressure calcula- 
tions for a single configuration of nuclei that we have 
obtain from DFT-MD with 57 atoms at 80 000 K. The 
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TABLE III: Table HI] continued. 





rs 


r(K) 


P (GPa) 


B/A/e (Ha) 


r~I / AT" / T T \ 

F/iVe (Ha) 


1. 


.86" 


1000 


13.55(3) 


-1.445347(12) 


-1.45574 


1, 


,86" 


3000 


21.50(8) 


-1.43837(4) 


-1.48081 


1, 


,86" 


5000 


28.04(9) 


-1.43196(4) 


-1.51097 


1, 


,86" 


10000 


41.40(10) 


-1.41740(6) 


-1.59608 
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zUUUU 


DO.lD^loj 


-l.ooDoo^y ) 


-i. / OOD 


1 
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9 91 /1Q 


I 


,86" 


finnnn 


167.7(2) 






1, 


,86" 


80000 


229.42(15) 


-1 04980fl2i 


-3.22147 


1 


.86^^ 


125000 


378(2) 


-0.743(3) 


-4.51072 


1 


.86*" 


250000 


918.6(1.3) 


0.2491(18) 


-8.65932 


1 


.86*" 


333333 


1340.7(1.3) 


0.9607(18) 


-11.7193 


1 


.86'' 


500000 


2214.5(1.7) 


2.336(2) 


-18.3468 


1 


.86'' 


571428 


2595(2) 


2.910(3) 


-21.3458 


1 


.86'' 


666667 


3104(3) 


3.661(4) 


-25.4472 


1 


.86'' 


800000 


3818(3) 


4.699(4) 


-31.3521 


1 


.86'' 


1x10*^ 


4876.7(1.5) 


6.212(2) 


-40.5306 


1 


.86'' 


2x10'' 


10128(3) 


13.559(3) 


-90.3193 


1 


.86'' 


4x10*^ 


20550(5) 


27.959(7) 


-199.554 


1 


.86'' 


8x10*^ 


41316(7) 


56.543(9) 


-437.845 


1, 


.86'' 


16x10'' 


82822 


113.69 


-955.409 


1, 


.86'' 


32x10** 


165822 


227.79 


-2068.45 


1 


.86'' 


64x10^ 


331768 


455.87 


-4450.87 


I 


ggd 


1 28 X 1 0^ 


663695 


911.93 


-9537.1 


1 

_L 


.ou 


oc;« V 1 


1 ^97^1 9 

-LOZr j O-LZr 


±0^4:. U-L 


9n'?'^s a 

-zjUOOO.O 


1 
1 


.oD 


c 1 w 1 fiS 
OIZ X lU 


ZD04DDO 


oD4o. 


/j Q1 Qf! /I 

-4oiyD.4 


1 
1 


.86 


1U24X lU 


ooUyoDD 


/ 296.30 


-914^^5.0 


1 


.86'' 


2048 XlO" 


10618754 


14592.66 


-193062 


1. 


,75" 


1000 


22.14(4) 


-1.443419(13) 


-1.45302 


1. 


,75" 


3000 


32.35(11) 


-1.43604(4) 


-1.47673 


1 

1, 


, 10 


IXAAA 

oUUU 


4U.DD(loj 


1 A onoo / K ^ 
-i.4zyoo(oj 


1 PiAPi'yi? 
-l.OUO/D 


1 

1. 


, 10 


iUUUU 


( .\y\j\V6) 


1 /1 1 /1 1 K 

-i.4i4io(Dj 


1 p;ooc7 
-l.OOOD/ 


1 

1. 


, ( 


zuuuu 






1 77/100 
-L.it 4ZZ 


I 


,75" 


4nnnn 






-9 IQ^SI 


1, 


,75" 


60000 


210.43(19) 


-1.18602(16) 


-2.66885 


1, 


,75" 


80000 


284.6(3) 


-1.05398(16) 


-3.18308 


1 


.75'' 


125000 


454.4(1.0) 


-0.7647(12) 


-4.44653 


1 


.75'' 


250000 


1098.0(1.1) 


0.2015(13) 


-8.50508 


1 


.75'' 


500000 


2639.1(1.6) 


2.2626(17) 


-17.9781 


1 


.75'' 


IxlO'' 


5831.7(2.0) 


6.143(2) 


-39.7109 


1 


.75'' 


2x10'' 


12139(2) 


13.503(3) 


-88.6228 


1 


.75'' 


4x10'' 


24656(3) 


27.915(4) 


-196.115 


1 


.75'' 


8x10'' 


49587(8) 


56.501(10) 


-430.926 


1, 


.75'' 


16x10'^ 


99420 


113.65 


-941.531 


I 


'jrd 

. 1 u 


32 X 10^ 




997 7(S 




1, 


.75" 


64x10'' 


398335 


455.85 


-4395.24 


1, 


.75" 


128x10*^ 


796789 


911.92 


-9425.79 


1, 


.75" 


256 xlO'' 


1593662 


1824.00 


-20115.4 


1, 


.75" 


512x10" 


3187384 


3648.11 


-42749 


1, 


.75" 


1024x10*^ 


6374809 


7296.30 


-90584.7 


1, 


.75" 


2048x10" 


12749648 


14592.65 


-191274 



nuclear coordinates are the — 1.86 snapshot are given 
in Tab. [X] in the appendix. For this fixed configuration 
of nuclei, we now compare the instantaneous pressure as 
a function of electronic temperature. The fact that the 
nuclei are now classical rather being represented by paths 
in PIMC has a negligible effect on the pressure for the 



TABLE IV: Table mil continued. 





T(K) 


P (GPa) 


E/Ne (Ha) 


F/Ne (Ha) 


1.5" 


1000 


73.92(10) 


-1.43440(3) 


-1.44135 


1.5" 
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1.5" 
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1.5'' 
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770(2) 


-0.7880(16) 
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1.5" 
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0.0907(12) 
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1.25" 
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331.6(3) 
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125000 


1565.4(1.5) 


-0.7817(4) 


-4.05965 


1.25'' 


250000 


3074(3) 


-0.0069(11) 


-7.65192 


1.25'' 


500000 


6999(4) 


1.8502(15) 


-15.9783 


1.25'' 


IxlO" 


15578(3) 


5.6796(11) 


-35.2469 


1.25'' 


2xl0" 


32897(7) 


13.107(3) 
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1.25'' 
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27.58(3) 
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113.381(13) 


-864.979 


1.25'' 
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227.54 


-1887.16 


1.25'' 
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1092767 


455.69 


-4087.74 


1.25'' 


128x10" 


2186203 


911.81 


-8810.14 


1.25'' 


256x10" 


4372877 


1823.92 


-18883.5 


1.25'' 


512x10" 
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3648.06 


-40286 


1.25'' 


1024x10" 


17492411 


7296.26 


-85659.8 


1.25'' 


2048x10" 


34984988 


14592.63 


-181419 



temperatures under consideration. In both methods, the 
instantaneous pressure is a well-defined quantity derived 
from the virial theorem. In DFT, one uses the diagonal 
elements of the stress tensor"^^, while is PIMC one derives 
the pressure directly from the kinetic, (/C), and potential, 
(V), energy, 

iPV = 2 (/C) + (V) , (1) 

where V is the volume of the simulation cell. DFT is pri- 
marily a groundstate method that use an exchange func- 
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TABLE V: Table ITVl continued. 



rs 


T(K) 


P (GPa) 


E/Ne (Ha) 


n /AT /T T \ 
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1" 
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1" 
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1* 
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1* 
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23334(5) 


3.8385(12) 


-24.9776 


1* 


IxlO'' 


29972('5~> 


5 3367^1 2~> 


-32 3609 


1* 
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1" 
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1" 
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-38656.7 
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1" 
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tional that was derived at r=0. It allows, however, to 
include electronic excitations using a thermal population 
of excited singe-particle states. On the other hand, PIMC 
is a finite temperature quantum simulation method that 
treats the electrons as interacting particles. The only 
approximation comes from the fermion nodes. 

Figure [4] compares the instantaneous pressures from 
both methods. At intermediate temperatures, there is a 
large interval where both methods agree. DFT pressures 
appear to be fairly accurate. For the level of accuracy 
needed for this study we could not detect any insuffi- 
ciency resulting from the groundstate exchange correla- 
tion functional nor from inaccurate thermal excitations 
resulting from an underestimated bandgap. However, the 
DFT eventually become prohibitively expensive at higher 
temperature. Some of the points at = 2.4 required up 
to 100 bands per atom, and that is one reason why we 
only used a single configuration. The other comes from 
path integrals. PIMC simulations with 123 atoms shown 
in Fig. [T] represent about the limit one can study with 
currently available computers. If one wants to repeat the 
calculations at lower temperature where the paths are 
longer, or using the more expensive variational nodes, 
one quickly exceeds existing limits in processing power. 

Figure [4] also reveals inaccuracies in the PIMC com- 
putation that are caused by approximations in the trial 
density matrix. PIMC with free-particle nodes predict 
pressures that are too high when the electrons settle into 
the ground state (T < 40 000 K for rs=1.86 as shown in 
Fig. This effect has already been reported for hy- 



FIG. 4: Comparison of the instantaneous pressure for a fixed 
configuration of nuclei derived from PIMC and DFT with 
thermally populated states. The upper panel included results 
from PIMC with variational nodes. The lower panel compares 
PIMC results with free-particle nodes for difi'erent densities. 
(The ideal gas contributions from the nuclei are not included.) 



drogenii. In the variational density matrix approach^^ 
one allows the trial density matrix to adjust to the po- 
sitions of the nuclei, which corrects most of the pressure 
error as can be seen in upper panel of Fig. [D However, 
the variational approach was derived to study finite tem- 
perature problems. It does not describe the electronic 
ground state as well as DFT. 

For the purpose of construction one EOS table for he- 
lium, we use DFT-MD results with electronic excitation 
up to 80 000 K for rs > 1.5 and results up to 125 000 K for 
rg — 1.25 and 1.0. For all higher T, we use PIMC simula- 
tions, which become more and more efficient at higher T 
because the length of the paths is inversely proportional 
to temperature. 

It should noted that the discussed validity range of dif- 
ferent trial density matrices depends very much on the 
material under consideration. The temperature where 
we switch from PIMC to DFT-MD reflects the degree 
of thermal electronic excitations as well as some depen- 
dence of the approximations made in both methods. The 
density dependence of the switching temperature would 
typically be estimated by comparing the temperature to 
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Temperature (K) 

FIG. 5: The relative excess pressure derived from PIMC 
(thick hnes with symbols) is compared with Debye-Hiickel 
plasma model (thin lines) for different values of rs given in 
the legend. The ideal pressure, Po, is derived Saha model of 
noninteracting helium species (see text). 

the Fermi energy of an ideal gas of electrons. However, 
to incorporate band structure effects of dense helium, we 
found it more appropriate to relate the switching temper- 
ature to the DFT band width. Band width and Fermi 
energy are identical in the systems of noninteracting par- 
ticles. For the purpose of this study, we found it appropri- 
ate to switch from PIMC to DFT-MD for temperatures 
corresponding to less than one third of the helium band 
width. 

We performed PIMC simulations up to 1.28 x 10^ K 
covering a large temperature interval of two orders of 
magnitude. The excess internal energies and pressures 
in Figs. [21 [31 and [5] are positive a lower temperatures 
reflecting electronic excitations but then change sign due 
to interactions of ions and free electrons. At very high 
temperature when helium is fully ionized, the system can 
be described by the Debye plasma model"^^. The Debye 
model is based on a self-consistent solution of the Poisson 
equation for a system of screened charges. The excess 
contribution to the free energy, internal energy, entropy 
per particle, and pressure are given by. 



Np 12 ' ATp 8 ' iVp 24 ' 24:V ' ^ ' 

i 

where ft = l/r^ is the inverse of Debye radius, rd, in a 
collection of Ni particles of charge Zi in volume V where 
A^p — J2i ^i- Figure [5] demonstrate very good agreement 
with the Debye model at high temperature. The De- 
bye model fails at lower temperatures where it predicts 
unphysically low pressures. Under these conditions the 
screening approximation fails because there are too few 



particles in the Debye sphere. The number of particles 
in the Debye sphere is proportional to, 

(r,/r,)3 ^ (Tr,)3/2 , (4) 

which means that the Debye model becomes increasingly 
accurate for high T and large r^. This is exactly what 
is observed in Fig. [5l For higher densities, PIMC and 
Debye predictions converge only at higher temperatures. 

The size of the Debye sphere increases with tempera- 
ture and will eventually exceed the size of any simulation 
cell. This occurs when the coupling the particles become 
very weak. With increasing temperature, the Coulomb 
energy decreases while the kinetic energy increase linearly 
with T . To determine the precise amount of the Coulomb 
energy becomes more and more difficult due to finite size 
effects. This is the reason why the PIMC simulations do 
not agree perfectly with the Debye model for the high- 
est temperature shown in the lower panel of Fig. [31 In 
conclusion, we use the Debye EOS for the highest data 
in our EOS table. The nonideal pressures reported in 
Fig. ini appear to be less sensitive to finite size errors than 
the internal energy because their volume dependence is 
relatively weak. 

Now we compare our first-principle EOS with chemi- 
cal free energy models that were developed before first- 
principles simulation data became available. Winisdoerf- 
fer and Chabrier'^'^ constructed a semianalytical model 
to study stellar interiors that covers a wide density range 
including metallization. Their EOS is only available in 
explicit form in a small temperature interval and a com- 
parison with DFT-MD simulation has already been re- 
ported^. That is why we focus on the free energy model 
derived by Saumon, Chabrier, and van Horn (SCvH)^. 
Together with their hydrogen model"*^, their EOS has 
been used numerous times to model giant planet interi- 
ors. 

Figure [HI compares the excess pressure from both EOSs 
for three different temperatures. At a very temperature 
of lO^K that are important for stellar interiors, we found 
fairly good agreement. The deviations between the SCvH 
model and PIMC simulations are only about 4%. 

At a intermediate temperature of 100 000 K, that ap- 
proximately represent the regime of shock wave experi- 
ments, the agreement is less favorable. One finds that the 
SCvH EOS reports pressures that are about 30% lower 
than those predicted by PIMC. This is partly due to the 
fact that the SCvH model follows the Debye model down 
to too low temperatures (Fig. [3]) and that an interpola- 
tion scheme between a low and a high-temperature model 
was used by SCvH that was not thermodynamically con- 
sistent (Fig. [131). 

The last panel in Fig. [51 is focused on the interiors of gi- 
ant planets with temperature of the order to 10 000 K. At 
low density both EOSs agree well, but above 1.5gcm~^ 
deviation begin to increase steadily. At conditions com- 
parable to Jupiter's interior, we find that the SCvH un- 
derestimates the pressure by 30%. In a hydrogen-helium 
mixture of solar composition, this translates into an er- 
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FIG. 6: Comparison of the relative excess pressure reported 
by first-principles simulations with the SCvH EOS model. 
The three temperatures shown here are relevant for stellar 
interiors, the comparison with shock wave experiments, and 
the interiors of giant planets. 

ror in the pressure of about 4%. This is the reason why 
even the helium EOS is important to estimate the size 
of Jupiter's core, which is expected to be only a small 
fraction of the Jupiter's total mass. 



IV. PAIR CORRELATION FUNCTIONS 

In this section, we study the structure of the fluid by 
analyzing correlations between different types of parti- 
cles. Given the large amount of simulation results, we 
focus our attention primarily on the temperature depen- 
dence and only report results for one density of = 1.86. 
The density dependence of the pair correlation functions, 
g{r), has been analyzed in Ref;~. 

Figure [7] shows how the nuclear pair correlation func- 




FIG. 7: The nuclear pair correlation functions for = 1.86. 
The lower panel shows the following temperatures: { 128, 64, 
32, 16, 8, 4, 2, 1, 0.5, and 0.125 } xlO*^ K from PIMC as weU 
as { 40, 20, 10, 5, 3, and 1 } xlO^ K from DFT-MD. 

tions changes over a temperature interval that spans 
seven orders of magnitude. At low temperature, the g{r) 
shows the oscillatory behavior that is typical for a hard- 
sphere fluid. The atomic interactions are governed by 
two tightly bound electrons that lead to a strong repul- 
sion at close range due to Pauli exclusion. As long as the 
density is not too high, this behavior is well-described by 
the Aziz pair potentiaUi. 

As temperature increases, two effects change the pair 
correlation function. The increase in kinetic energy leads 
to stronger collisions and atoms approach each other 
more. This is regard, helium is not exactly a hard-sphere 
fluid because the Aziz pair potential is not perfectly hard. 
The increase in temperature also damps of the oscillation 
in the g{r). 

At 80 000 K, one finds perfect agreement between 
PIMC and DFT-MD (upper panel in Fig. [T]). As the 
temperature is increased further, changes in nuclear g(r) 
function are dominated by electronic excitations and the 
ionization of atoms. One finds that the strong repul- 
sion at low temperature disappears gradually. As one 
approaches the Debye-Hiickel limit, the fiuid behaves like 
correlated system of screened Coulomb charges. 

The peak in the electron-nucleus pair correlation func- 
tions in Fig. [5] illustrates that the electrons are bound 
to the nuclei. At 40 000 K and below, the peak height 
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FIG. 8: The electron-nucleus pair correlation functions, g{r), 
from PIMC for rs = 1.86. Starting with the highest peak, the 
following temperatures are plotted: { 0.04, 0.08, 0.125, 0.25, 
0.333, 0.5, 0.8, 1, 2, 4, and 128 } x lO*' K. We plot r * g{r) on 
the ordinate so that the peak at small r illustrates the fraction 
of electrons in bound states. The decrease in peak height 
with increasing temperature demonstrates thermal excitation 
of electrons, which eventually leads the ionization of atoms. 




FIG. 9: The electron-electron pair correlation functions for 
electrons with parallel spins calculated with PIMC for Vs = 
1.86. Starting from the left, the following temperatures are 
plotted: { 32, 16, 8, 4, 2, 1, 0.25, 0.125, 0.08, 0.0625, 0.04 } 
xlO^K. With increasing temperature, correlation effects are 
reduced and the exchange-correlation hole disappears. 



is maximal. At higher temperature, electrons get ex- 
cited and eventually atoms get ionized. The peak height 
is consequently reduced until, at very high temperature, 
the motion of electrons and nuclei appears to be uncor- 
related. 

The correlation of electrons with parallel spins is de- 
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FIG. 10: The electron-electron pair correlation functions for 
electrons with opposite spins calculated with PIMC for = 
1.86. Starting from the top, the following temperatures are 
plotted: { 0.04, 0.0625, 0.08, 0.125, 128, and 1 } xlO'^K. The 
smallest values are observed for 10® K. 

termined by Pauli exclusion and Coulomb repulsion but 
is also influenced by the motion of the nuclei at low tem- 
perature. Combination of all these effects causes the mo- 
tion of same-spin electrons to be negatively correlated 
at small distances. This is typically referred to as the 
exchange-correlation hole. At high temperatures, kinetic 
effect reduces the size of this hole but g{r) always goes 
to zero for small r due to Pauli exclusion. 

Despite the Coulomb repulsion, the electrons with op- 
posite spins are positively correlated at low temperature, 
because two electrons with opposite spin are bound in a 
helium atom. With increasing temperature, the peak in 
Fig. [TO] reduces in height because more and more elec- 
trons get ionized. At 10®, one finds the lowest values for 
g{r — > 0) because the electrons are anti-correlated due to 
the Coulomb repulsion. If the temperature is increased 
further kinetic effects dominate over the Coulomb repul- 
sion and g(r — > 0) again increases and will eventually 
approaches 1 at high temperature. 

V. ENTROPY CALCULATIONS 

Convection in the interior of planets determines that 
the temperature-pressure profile is adiabatic. In con- 
sequence, the planetary interiors is fully determined by 
the conditions on the surface and the EOS. This makes 
the calculation of adiabats important. However, neither 
Monte Carlo nor molecular dynamics methods can di- 
rectly compute entropies because both techniques save 
orders of magnitude in computer time by generating a 
only representative sample of configurations. Without 
this gain many-body simulations would be impossible 
and in consequence, entropies that are measure of the 
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FIG. 11: Comparison in temperature-density space of adi- 
abats from first-principles simulations (this work) and the 
SCvH EOS model. 



total available phase space are not accessible directly. 

Typically one derives the entropy by thermodynamic 
integration from a know reference state. However, for the 
planetary interiors, the absolute value of the entropy is 
not important as long as one is able to construct (T, P) 
curve of constant entropy. This can be done using the 
pressure and the internal energy from first-principles sim- 
ulations at different (T, V) conditions. Using Maxwell's 
relations, one finds. 
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(5) 



By solving this ordinary differential equation, (V,T)- 
adiabats can be constructed as along as a sufficiently 
dense mesh of high-quality EOS points are available to 
make the required interpolation and differentiation of E 
and P with respect to temperature satisfactorily accu- 
rate. 
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FIG. 12: Comparison in temperature-pressure space of the 
adiabats shown in Fig. 1111 



One drawback of formula (O is that it is not necessarily 
thermodynamically consistent if pressures and internal 
energies are interpolated separately. This is the primary 
reason why we developed the following method to fit the 
free energy instead. Pressure and internal energy are 
related to the free energy, F{V,T), by 



P = - 



dF 
dV 



and 



E = F-T 



dF 
df 



(6) 



Different EOS fits for fluids have been proposed in the 
literature^i^. Thermodynamic consistency was not a 
priority in each case. Both papers relied on specific func- 
tional forms that were carefully adjusted to the material 
under consideration. Although such a fit of specific form 
could probably have also been constructed for the pre- 
sented helium EOS data, we wanted to have an approach 
that is not just applicable to one material. Therefore, we 
decided to represent the free energy as a bi-cubic spline 
function with temperature and density as parameters. 
This spline function can accurately represent our helium 
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EOS data and can easily be adapted to fit other mate- 
rials. Cubic splines are twice continuously difFerentiable, 
which means the derived pressures and energies are once 
continuously difFerentiable with respect to V and T. This 
is sufficient for this study. If additional thermodynamic 
functions that require higher order derivatives of the free 
energy, such as sound speeds, need to be fit also then 
higher order splines can accommodate that. 

We start the free energy interpolation by constructing 
as series of one-dimensional splines functions Fn{T) for 
different densities. The choice of knots Tj is arbitrary. 
Their location should be correlated with the complex- 
ity of the EOS as well as the distribution of EOS data 
points. In our helium example, we used a logarithmic 
grid in temperature with about half as many knots as 
data points. The set of free energy values on the knots, 
F{Ti), represent the majority of the set of fit parameters. 
In addition, one may also include the first derivatives of 
the splines ^^\v the lowest and highest temperature, 
which represent the entropy. Alternatively, one could de- 
rive those derivatives by other means and then keep them 
fixed during the fitting procedure. 

To compute the free energy at a specific density, n*, 
and temperature, T*, we first evaluate all sphnes F„(T*) 
and then construct a secondary spline at constant tem- 
perature as function of density. Ft- (n) . Its first derivate 
is related to the pressure. Again, the derivative at the 
interval boundaries can either be fixed or adjusted during 
the fitting procedure. We adjust them by introducing an 
additional spline ^\t{T) at the lowest and highest den- 
sities, which then get adjusted in the fitting procedure. 

We begin the fitting procedure with an initial guess 
for the free energy function derived from Eq. Then 
we employ conjugate gradient methods^^ to optimize the 
whole set of fitting parameters. Minimizing the sum of 
the squared relative deviations in pressure and internal 
energy has been found to work best. (Just for the deriva- 
tion of the relative deviation in energy, the zero of energy 
has been shifted to value of the isolated helium atom.) 

All fits tend to introduce wiggles if too many free pa- 
rameters are included. We control wiggles by adjusting 
the number of knots in density and temperature but we 
also introduced penalty in the form. 




to favor fits with smaller \d^P/dn^\. Finally we changed 
the density argument in the spline interpolation from 
Ft{p) to FT(log(p)). This improves the fit in the high 
temperature limit where the free energy is dominated by 
the ideal gas term that has logarithmic dependence on 
density. 

The presented free energy fit is thermodynamically 
consistent by construction. It allowed us to accurately 
represent the entire data set of P and E values. Without 
additional information, the free energy can be determined 
up to a term TAS, which is sufficient to compute adia- 
bats. To determine the absolute value of the entropy, one 




10"' 10° 10' 10' 10^ 10' 10^ 
Pressure (GPa) 

FIG. 13: Comparison of Jupiter's isentrope with helium shock 
Hugoniot curves for different precompression ratios. The la- 
bels specify the precompression ratio relative to the density 
at ambient pressures (po ~ 0.1235 g cm~'^, Ref.— ). The sym- 
bols approximately represent recent experiments^^. The in- 
side of the dashed box indicates conditions, for which the 
SCvH EOS^ was interpolated and is not thermodynamically 
consistent. 

needs an anchor point, for which the entropy was derived 
by different means. 

Figure [TT] compares different adiabats derived from our 
first-principles EOS with predictions from the SCvH EOS 
model. Beginning from a joint starting point of = 2.4 
and a selection of seven different temperatures of 3000, 
5000, 10 000, 50 000, 100 000, 500 000, and 10*^ K, we con- 
structed the adiabats for both models for the density in- 
terval under consideration. The upper panel of Fig. [11] 
demonstrates good agreement between both methods a 
low densities up to about 1 g cm~^. For higher densi- 
ties, one finds deviations of up to 20% in the predicted 
temperatures on the adiabat. A higher temperature, the 
agreement get substantially better, which is illustrated 
in the lower panel of Fig. [TTJ The observed deviations 
are similar to pressure differences shown in Fig. [SI 

For applications in the field of planetary science, we 
also show the adiabats in {T,P) space in Fig. [121 The de- 
viations are comparable in magnitude but appear smaller 
on a logarithmic scale. 



VI. SHOCK WAVE EXPERIMENTS 

Dynamic shock compression experiments are the pre- 
ferred laboratory experiments to probe the properties 
materials at high pressure and temperature. Using big- 
ger and more powerful lasers, great progress has been 
made in reaching higher and higher pressures. Under 
shock compression, the initial state of material character- 
ized by internal energy, pressure, and volume {Eq, Pq, Vq) 
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FIG. 14: Hugoniot curves for different precompression ratios 
from Fig. [13] plotted as function of shock compression. 



changes to the final state described by {E,P,V). The 
conservation of mass, momentum, and energy yields the 
Hugoniot condition*^^ , 

H={E-Eo) + ^{P + Po){V-Vq) = 0. (8) 

Different shock velocities lead to a collection of final 
states that are described by a Hugoniot curve. Using 
Eq. [3 this curve can easily calculated for a given EOS 
where one most often may assume Pa <C P. Vo = 32.4 
cm'^/mol {pq = 0.1235 gcm~^) is taken from experi- 
ment^^. For Ef), one takes the energy of an isolated 
helium atom, which must be calculated consistently with 
the final internal energy, E. An initial static precompres- 
sion that changes Vq will also affect Eq and Pq but the 
corrections are negligible as long as the amount of initial 
compression work as small compared to the energy that 
is deposited dynamically. Assuming cIEq — dPo = 0, the 
total differential of H reads. 



dH = dE + -dV - -dVo 
2 2 



^{V-VQ)dP (9) 



The point of maximum compression, T^max = Vq/V, 
along the Hugoniot curve can be derived by setting, 
dH = dV = dVa = 0. The resulting condition can 
be expressed in terms of the Griineisen parameter, 7 = 

^iilv=2/('7max-l). 

Now we will determine how the maximum compres- 
sion ratio, ?7max, changes if the sample is precompressed 
statically. Keeping the final shock pressure constant, the 
compression ratio changes as function of the initial sam- 
ple volume, Vq, 



drj 



1 Vq dV 

V^Y^Wq 



(10) 
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FIG. 15: Pressure-density plot of shock Hugoniot curves for 
different precompression ratios. The upper panel shows the 
shock Hugoniot curves from Fig. 1141 recent experimental re- 
sults (symbols), and the range of the lower plot. Below we 
plot the Hugoniot curves for precompression ratios (see la- 
bels) that approximately match the experimental conditions 
(symbols, see Ref.— ). The error bar on the upper left solid 
line represents the uncertainty in the calculations. 



Setting diJ = dP = in Eq. ^ one finds. 



dV 



2 dE 
PdV 



Inserting this result into Eq. [TO] yields, 



V 



2(7 - 5) 
7(2 + <5) 



(11) 



(12) 



Since the parameters 7 and 5 are both positive, the re- 
> 0, is equivalent to the relation. (5 < 7, 



lation, 



dVa 



which is again equivalent to, 

P_ ' 
P 'dp 



(13) 
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If this condition is fulfilled for a particular EOS then 
the maximum shock compression ratio will decrease if 
the sample is precompressed statically, which reduces Vq. 
We have computed the isoenergetic compressibility for 
our first-principles EOSs for helium and hydrogen and 
verified that this condition is satisfied for both materials 
(Fig. [T4|) . It is also fulfilled for an ideal plasma model be- 
cause the maximum compression ratio is determined by 
the balance of excitations of internal degrees of freedom 
and interaction effectsi^^. Although all interactions are 
neglected, an ideal model correctly represents the fact 
the excited states are suppressed at high density because 
of the reduced entropy. The diminished importance of 
excitations reduces the maximum compression ratio to 
values closer to 4, which is the expected result for non- 
interacting systems without internal degrees of freedom. 

Recent laser shock wave experiments^^ reached pres- 
sures of 2 megabars in fluid helium for the first time. 
The sample was precompressed statically in a modified 
diamond anvil cell before the shock was launched. The 
static precompression is an important development that 
enables one to reach higher densities and still allows one 
to direct determine the EOS. Reaching higher densities is 
important for planetary interiors because shock Hugoniot 
curves rise faster than adiabats in a P-T diagram shown 
in Fig. [131 As a result, a large part of Jupiter's adia- 
bat remains inaccessible unless one increases the start- 
ing density by precompression. The precompression and 
relation of planetary interiors was studies theoretically 
in Ref.''^. It was demonstrated that precompression of 
up to 60 GPa would be needed to characterize 50% of 
Jupiter's envelope. The challenge is here to reach high 
enough densities because a single shock wave compresses 
the material only 5.25-fold or less (Fig. [T4|) . 

The measurements of J. Eggert et al}& confirmed two 
of our theoretical predictionsi^. They showed that he- 
lium has a shock compression ratio substantially larger 
than 4 due the electronic excitations and that the com- 
pression ratio would decrease with increasing precom- 
pression (Fig. ll4p . However, the measurements appeared 
to be in better agreement with the SCvH EOS model 
than PIMC simulations^^. The SCvH model predicts a 
maximum compression ratio of 6.5 to occur around 300 
GPa. A different chemical model based on an expansion 
of the activity*^ predicts maximum compression ratios 
between 5.6 and 6.2 to occur at about 100 GPa. 

Figure [151 shows a detailed comparison between exper- 
iments and our first-principles simulations. The shock 
measurements without precompression indeed show a 



higher compression than predicted from first principles. 
The deviations are outside the experimental error bars. 
However, this discrepancy goes away with increasing pre- 
compression. The shocks with 3.4-fold precompression 
are in good agreement with first-principles predictions. 
We have no explanation for this trend at present. 

The reason why the SCvH EOS yields larger compres- 
sion ratios can be understood by looking at the pressure 
that this model predicts. Using our first-principles EOS, 
we derived the shock temperatures that correspond the 
reported P-p measurements. In the resulting tempera- 
ture range of 24 000 - 63 000 K, the SCvH EOS signif- 
icantly underestimates the pressure (see Fig. [6]), which 
leads to higher predicted compressions (Eq. [8]). Further- 
more, all measurements fall in the region where the SCvH 
model relied on interpolation (Fig. I13p and is not expect 
be as reliable. 



VII. CONCLUSIONS 

This paper combined path integral Monte Carlo and 
density functional molecular dynamics simulation to de- 
rive one coherent equation of state for fluid helium at 
high pressure and temperature. Helium is a compara- 
tively simple material since it does not form chemical 
bonds nor has core electrons, but the our approach of 
combining two simulation techniques can be generalized 
to study more complex materials at extreme conditions. 
Certainly the presented approach to fit the free energy 
and to derive adiabats works for any set of EOS data 
points derived from first-principles simulations. 

For the future, one might consider replacing DFT-MD 
with coupled ion-electron Monte Carlc4^. However this 
is strictly a groundstate method and one would still need 
to find a way to include electronic excitations. 
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APPENDIX A: FREE ENERGY SPLINE 
INTERPOLATION 



We constructed the following 2D spline interpolation 
of the free energy in order to reproduce the internal en- 
ergy and pressures from Tab. [D We use atomic units 
of Hartrees and Bohr radii. For each density of = 
{2.4, 2.0, 1.6, 1.2, 0.8}, we construct a cubic spline i^„(T). 
Table EJ lists 16 knot points {T„F{T,)) for each den- 
sity. In addition, the first derivate ^ are specified at 
the lowest and highest temperatures. This is sufficient to 
construct a cubic spline function _F(T)'*^. 

In a similar fashion, we derive a spline function that 
contains that free energy derivative with respect to den- 
sity, ^(T), at the lowest and highest densities, = 2.4 
and 0.8 respectively, n is the density of the electrons, 
n = Ne/V. Those knot points as well as the T deriva- 
tives are included in Tab. IVII also. 

In order to obtain the free energy for a particular 
density and temperature, (n*,T*), we proceed as fol- 
lows. First we evaluate the spline functions F{T*) and 
^(T*) at temperature T*. Using these five knots points 
and density derivatives, we construct a spline function, 
F(log(n)). We use log(n) as argument because it better 
represents the high-temperature limit of weak interac- 
tions. Note that the constructed splines for the density 
derivate contain ^ and not jh^^-t- Then i^(log(n)) is 

on alog(n) ^ 

evaluated at the density of interest, n* . Finally we add 
the term, -TAS = -13.7902836 Ha*T, which brings 
the entropy in agreement with our Debye-Hiickel refer- 
ence point at high temperature for rs=1.86. This pro- 
cedure yields the free energy F(n*,T*) in Hartrees per 
electron. Other thermodynamic variables including pres- 
sure, internal energy, entropy, and Gibbs free energy can 
be obtained by differentiation. 
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TABLE VI: Knot points for free energy spline interpolation 
r(a.u.) f(r, ^ 2.4, T) /(r, = 2.0, T) /(r, = 1.6, T) /(r, ^ 1.2, T) /(r, ^ 0.8, T) 



0.001583407607 
0.004369348882 
0.01205704051 
0.03327091284 
0.09180973061 
0.2533452171 
0.6990958211 
1.929126481 
5.323346054 
14.6895569 
40.53523473 
111.8553314 
308.660237 
851.7353686 
2350.329104 
6485.637557 
/'(r.,Ti) 
f'{rs,TN) 



-1.433567121 
-1.406243368 
-1.338444742 
-1.169814789 
-0.7578793569 
0.1213432393 
1.448161918 
1.650153009 
-6.480107758 
-50.37972679 
-231.4313924 
-893.7314032 
-3172.519712 
-10693.37579 
-34893.6373 
-111050.5942 
10.43646526 
-19.3620206 



-1.431135811 
-1.402248169 
-1.330719999 
-1.154154406 
-0.7214892496 
0.2279596473 
1.849742465 
3.082886945 
-2.223161966 
-38.41515844 
-198.2806566 
-802.1449141 
-2918.919455 
-9996.590015 
-32971.41578 
-105746.9466 
10.93007841 
-18.54727965 



-1.422237604 
-1.390358674 
-1.313137294 
-1.125146478 
-0.6653645937 
0.3697481593 
2.331720062 
4.795803813 
2.964697691 
-23.79211621 
-157.6765535 
-689.7861073 
-2608.696846 
-9137.681211 
-30600.55572 
-99203.62441 
12.42958105 
-17.53417624 



-1.377317056 
-1.34121764 
-1.255091313 
-1.047701304 
-0.5501514183 
0.5961317377 
2.967152398 
6.956651264 
9.579180258 
-4.992619324 
-105.3657984 
-544.9061891 
-2207.77895 
-8032.176108 
-27552.36594 
-90788.1102 
14.04097874 
-16.23883073 



-1.089747214 
-1.049220122 
■0.9491949603 
■0.7135001905 
■0.1583196636 
1.145222958 
4.011622282 
9.947352284 
18.74636894 
21.39894481 
-31.61456362 
-340.7935545 
-1644.63782 
-6477.239292 
-23259.81241 
-78945.99564 
13.75171132 
-14.41224513 



fir,, y^'^ 



r(a.u.) 



: 2.4, T) 



0.001583407607 

0.004369348882 
0.01205704051 
0.03327091284 
0.09180973061 
0.2533452171 
0.6990958211 
1.929126481 
5.323346054 
14.6895569 
40.53523473 
111.8553314 
308.660237 
851.7353686 
2350.329104 
6485.637557 

/'(r-.,ri) 



0.1605721538 

0.2923651353 
0.6458549976 
1.463747273 
3.574255133 
10.98479122 
43.27443042 

153.211445 
450.4159975 
1263.720981 
3502.359609 
9679.619949 
26719.83852 

73736.1465 
203474.5198 
561480.8894 
47.35562925 
86.57538487 



0.9460950728 
0.9625433841 
0.9782855175 
1.029890765 
1.142591093 
1.409424435 
2.089943658 
5.309130852 
15.9843674 
46.33079487 
129.8984992 
359.8256334 
992.5003633 
2741.461026 
7568.332998 
20879.15481 
6.453623893 
3.218499127 
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TABLE VII; Reduced coordinates of the DFT-MD configu- 
ration with 57 atoms that was used to report the = 1.86 
results for the instantaneous pressure in Fig. |4] The cell size 
is L = 14.5382 a.u. 





x/L 




y/L 




z/L 




x/L 




y/L 




z/L 





749029 





334272 





723992, 





359050 





631169 





090795 





636183 





917961 





531890, 





500277 





715818 





420142 





509121 





642554 





328933, 





192192 





222632 





042651 





273631 





845722 





363632, 





070837 





830223 





693497 





053785 





837401 





054990, 





138489 





091713 





097622 





250609 





517490 





740851, 





953625 





430789 





067921 





107008 





407958 





463387, 





023708 





960709 





487179 





988548 





830572 





241931, 





811738 





062550 





902069 





244399 





482412 





399190, 





693258 





647174 





360832 





924284 





678572 





470508, 





181701 





886709 





333868 





780287 





033015 





620919, 





859185 





932541 





252564 





774645 





083064 





349744, 





903457 





888628 





124621 





293881 





081041 





053630, 





220134 





760599 





688370 





493690 





930407 





343378, 





585411 





439278 





167284 





648043 





965342 





702852, 





219455 





957094 





895428 





504966 





639074 





084498, 





906610 





508304 





938057 





716468 





854022 





986517, 





385839 





307391 





681601 





099368 





291429 





740170, 





475139 





160612 





598743 





252564 





696499 





576596, 





788211 





564812 





486616 





613177 





259980 





238984, 





296858 





344416 





229757 





526564 





816547 





598836, 





429733 





712523 





742929 





507514 





904602 





268688, 





685066 





562001 





926251 





614731 





263859 





402947, 





432246 





210193 





939664 





115992 





498747 





676389, 





424152 





141821 





676522 





778767 





981750 





935757, 





208696 





768371 





292528 





334815 





183086 





275601, 





487257 





590889 





227333 





975542 





456665 





257836, 





577884 





835181 





876629 





737370 





699890 





544111, 





177496 





781162 





853225 





558513 





066648 





194491 















